Build gene co-expression network (GCN) from time-course gene expression data

---

Last updated on 2025-09-04.

---

library(dplyr)
library(dbplyr)
library(ggplot2)
for (i in list.files(here::here("R"), full.names = TRUE)) {
  source(i)
}

# SAMPLE NAME
## specify sample name
sample.names <- c(
  # dmel
  "dmel-head",
  "dmel-body",
  "dmel-heart"
  # # mmus
  # "mmus-brain_stem", 
  # "mmus-cerebellum", 
  # "mmus-hypothalamus", 
  # "mmus-heart", 
  # # panu
  # "panu-brain_stem",
  # "panu-cerebellum",
  # "panu-hypothalamus",
  # "panu-heart",
  # "panu-scn"
)
# sample.cycles <- c("LD", "DD")

## SPECIFY THE DATASET TO BUILD GCN WITH
which.sample <- "dmel-heart"

writeLines(
  glue::glue("Sample: {which.sample}")
)

Sample: dmel-heart

Structure of input data:

ONE-SHOT (make_modules)

Make and plot modules

mods <- purrr::map(
  sample.names,
  ~ timecourseRnaseq::make_modules(
      data = tidydata.db[[.x]],
      log2 = TRUE,
      id_column = "gene_name",
      min_expression = NULL,  # automatically estimated
      min_timepoints = NULL,  # automatically estimated
      method = "wgcna",
      qc = TRUE,
      sim_method = "kendall",
      soft_power = 15,      # automatically estimated
      min_module_size = 50,
      max_modules = 16,
      merge_cutoff_similarity = 0.9,
      plot_network_min_edge = 0.5,
      plot_network = FALSE,
      tidy_modules = TRUE
    )
) |>
  setNames(
    sample.names
  )

---------------------------------------------------
1. Log2-transform and subset 
---------------------------------------------------
Applying log2-transformation...Done.
Estimated min_expression = 4 
Estimated min_timepoints = 6 
Subsetting data...Done.
[ NOTE ]: After subsetting, 3100 of 8212 rows remain.

 Flagging genes and samples with too many missing values...
  ..step 1
All okay!Visualizing the log-transformed data
---------------------------------------------------
2. Calculate similarity of expression 
---------------------------------------------------
Running gene-gene similarity...Done.


---------------------------------------------------
3. Create adjacency matrix 
---------------------------------------------------
Performing network topology analysis to pick
  soft-thresholding power...
   Power SFT.R.sq    slope truncated.R.sq mean.k. median.k. max.k.
1      1 0.967000  2.94000          0.958    1610    1690.0   2010
2      2 0.956000  1.26000          0.944    1090    1140.0   1580
3      3 0.918000  0.66200          0.895     827     854.0   1320
4      4 0.876000  0.33900          0.843     663     670.0   1160
5      5 0.664000  0.13700          0.578     551     551.0   1030
6      6 0.000259 -0.00156         -0.280     469     467.0    936
7      7 0.611000 -0.10600          0.501     408     401.0    858
8      8 0.757000 -0.19700          0.690     359     350.0    794
9      9 0.864000 -0.26400          0.829     320     309.0    739
10    10 0.869000 -0.31500          0.839     288     274.0    692
11    12 0.877000 -0.40500          0.854     238     221.0    614
12    15 0.873000 -0.49700          0.865     187     164.0    527
13    18 0.874000 -0.57000          0.884     152     127.0    461
14    21 0.875000 -0.62700          0.890     126      99.7    410

Plotting resutls from the network topology analysis...

Done.
[ NOTE, FIGURE ]: Red horizontal line indicates a signed R^2 of 0.9

Setting soft-thresholding power to: 15.
Power-transforming the gene-gene similarity matrix...Done.

---------------------------------------------------
4. Convert into topological overlap matrix (dissTOM) 
---------------------------------------------------
Creating dissTOM...Done.
Performing hierarchical clustering on dissTOM...Done.

---------------------------------------------------
5. Identify modules (clusters) 
---------------------------------------------------
Merging modules that have a correlation ≥ 0.9 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
         blue          cyan     darkgreen      darkgrey    darkorange 
          843           173           342           174            66 
      darkred darkturquoise         green    lightgreen   lightyellow 
           82            78           257            90           347 
      magenta     turquoise        yellow 
          251           234           163 

Cutoff used: 0.9 
Number of modules identified: 13 

Calculating module-module similarity based
  on module-eigengene-expression...Done.
Tidying module names...Done.
Plotting adjacency matrix for module-module similarity...

---------------------------------------------------
6. Tidy modules (clusters) 
---------------------------------------------------
---------------------------------------------------
1. Log2-transform and subset 
---------------------------------------------------
Applying log2-transformation...Done.
Estimated min_expression = 4 
Estimated min_timepoints = 6 
Subsetting data...Done.
[ NOTE ]: After subsetting, 3581 of 8210 rows remain.

 Flagging genes and samples with too many missing values...
  ..step 1
All okay!Visualizing the log-transformed data
---------------------------------------------------
2. Calculate similarity of expression 
---------------------------------------------------
Running gene-gene similarity...Done.


---------------------------------------------------
3. Create adjacency matrix 
---------------------------------------------------
Performing network topology analysis to pick
  soft-thresholding power...
   Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1  0.76400  6.240          0.970  1300.0   1310.00 1610.0
2      2  0.59600  2.700          0.974   667.0    676.00  971.0
3      3  0.31200  1.220          0.961   403.0    403.00  672.0
4      4  0.05370  0.378          0.932   267.0    263.00  500.0
5      5  0.00494 -0.103          0.927   189.0    183.00  391.0
6      6  0.07690 -0.393          0.944   139.0    133.00  316.0
7      7  0.19200 -0.622          0.957   106.0     99.50  261.0
8      8  0.31400 -0.811          0.966    83.5     76.60  220.0
9      9  0.41400 -0.955          0.972    67.0     60.20  189.0
10    10  0.48500 -1.060          0.977    54.7     48.30  164.0
11    12  0.60600 -1.240          0.986    38.1     32.80  127.0
12    15  0.71900 -1.410          0.991    24.1     19.80   91.6
13    18  0.78600 -1.510          0.997    16.3     12.80   69.2
14    21  0.81200 -1.600          0.993    11.6      8.79   54.2

Plotting resutls from the network topology analysis...

Done.
[ NOTE, FIGURE ]: Red horizontal line indicates a signed R^2 of 0.9

Setting soft-thresholding power to: 15.
Power-transforming the gene-gene similarity matrix...Done.

---------------------------------------------------
4. Convert into topological overlap matrix (dissTOM) 
---------------------------------------------------
Creating dissTOM...Done.
Performing hierarchical clustering on dissTOM...Done.

---------------------------------------------------
5. Identify modules (clusters) 
---------------------------------------------------
Merging modules that have a correlation ≥ 0.9 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
        black          blue         brown          cyan     darkgreen 
          149           190           181           118            90 
     darkgrey    darkorange       darkred darkturquoise         green 
           79            65            91            88           165 
  greenyellow        grey60     lightcyan    lightgreen   lightyellow 
          133           105           108           104           104 
      magenta  midnightblue        orange paleturquoise          pink 
          144           112            71            50           148 
       purple           red     royalblue   saddlebrown        salmon 
          134           163            96            60           118 
      skyblue     steelblue           tan     turquoise         white 
           61            52           121           238            63 
       yellow 
          180 
Merging modules that have a correlation ≥ 0.85 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
        black          blue         brown          cyan     darkgreen 
          149           190           181           118           223 
     darkgrey    darkorange       darkred darkturquoise         green 
          200            65            91            88           165 
       grey60     lightcyan    lightgreen   lightyellow       magenta 
          105           108           167           104           144 
 midnightblue        orange paleturquoise          pink        purple 
          112            71            50           148           134 
          red     royalblue   saddlebrown        salmon       skyblue 
          163            96            60           118            61 
    steelblue     turquoise        yellow 
           52           238           180 
Merging modules that have a correlation ≥ 0.8 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
        black          blue         brown          cyan     darkgreen 
          240           190           269           118           461 
     darkgrey    darkorange         green        grey60     lightcyan 
          200            65           165           105           108 
   lightgreen   lightyellow       magenta  midnightblue        orange 
          167           104           196           112            71 
paleturquoise          pink        purple           red     royalblue 
          168           148           134           163            96 
  saddlebrown       skyblue        yellow 
           60            61           180 
Merging modules that have a correlation ≥ 0.75 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
        black         brown          cyan     darkgreen    darkorange 
          240           269           118           857            65 
       grey60     lightcyan    lightgreen   lightyellow        orange 
          105           108           634           104            71 
paleturquoise          pink        purple           red     royalblue 
          168           148           134           163           157 
  saddlebrown        yellow 
           60           180 
Merging modules that have a correlation ≥ 0.7 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
      black        blue       brown        cyan   darkgreen  darkorange 
        240         634         403         118         857          65 
     grey60   lightcyan lightyellow      orange        pink         red 
        105         276         104          71         148         163 
  royalblue saddlebrown      yellow 
        157          60         180 

Cutoff used: 0.7 
Number of modules identified: 15 

Calculating module-module similarity based
  on module-eigengene-expression...Done.
Tidying module names...Done.
Plotting adjacency matrix for module-module similarity...

---------------------------------------------------
6. Tidy modules (clusters) 
---------------------------------------------------
---------------------------------------------------
1. Log2-transform and subset 
---------------------------------------------------
Applying log2-transformation...Done.
Estimated min_expression = 3 
Estimated min_timepoints = 3 
Subsetting data...Done.
[ NOTE ]: After subsetting, 3461 of 8258 rows remain.

 Flagging genes and samples with too many missing values...
  ..step 1
All okay!Visualizing the log-transformed data
---------------------------------------------------
2. Calculate similarity of expression 
---------------------------------------------------
Running gene-gene similarity...Done.


---------------------------------------------------
3. Create adjacency matrix 
---------------------------------------------------
Performing network topology analysis to pick
  soft-thresholding power...
   Power SFT.R.sq    slope truncated.R.sq mean.k. median.k. max.k.
1      1 7.94e-01  3.84000        0.82900    1820      1840   2130
2      2 6.57e-01  1.55000        0.65800    1260      1260   1610
3      3 5.03e-01  0.82900        0.52500     964       957   1320
4      4 2.77e-01  0.45400        0.29500     783       771   1130
5      5 1.27e-01  0.26000        0.13800     661       644    986
6      6 2.56e-02  0.10300        0.00816     572       552    880
7      7 4.07e-05  0.00404       -0.01370     505       483    796
8      8 2.55e-02 -0.09490        0.04560     451       429    728
9      9 7.55e-02 -0.16300        0.10800     409       385    671
10    10 1.35e-01 -0.22000        0.16200     373       349    622
11    12 2.29e-01 -0.29500        0.27700     318       294    545
12    15 3.18e-01 -0.37100        0.36900     261       238    460
13    18 3.82e-01 -0.42500        0.44600     221       200    399
14    21 4.16e-01 -0.45400        0.48600     192       173    353

Plotting resutls from the network topology analysis...

Done.
[ NOTE, FIGURE ]: Red horizontal line indicates a signed R^2 of 0.9

Setting soft-thresholding power to: 15.
Power-transforming the gene-gene similarity matrix...Done.

---------------------------------------------------
4. Convert into topological overlap matrix (dissTOM) 
---------------------------------------------------
Creating dissTOM...Done.
Performing hierarchical clustering on dissTOM...Done.

---------------------------------------------------
5. Identify modules (clusters) 
---------------------------------------------------
Merging modules that have a correlation ≥ 0.9 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
        black          blue         brown          cyan     darkgreen 
          178           252           238           126            66 
      darkred darkturquoise         green   greenyellow        grey60 
           69            53           220           141            93 
    lightcyan    lightgreen   lightyellow       magenta  midnightblue 
           93            87            80           172            99 
         pink        purple           red     royalblue        salmon 
          175           170           214            71           131 
          tan     turquoise        yellow 
          137           364           232 
Merging modules that have a correlation ≥ 0.85 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
        black          blue         brown          cyan     darkgreen 
          542           424           413           126            66 
      darkred darkturquoise         green   greenyellow        grey60 
           69            53           220           141            93 
    lightcyan    lightgreen   lightyellow  midnightblue        purple 
           93           224            80            99           301 
          red     royalblue        yellow 
          214            71           232 
Merging modules that have a correlation ≥ 0.8 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
        black          blue         brown          cyan     darkgreen 
          542           424           413           126           159 
      darkred darkturquoise         green   greenyellow     lightcyan 
          149            53           434           373           164 
   lightgreen  midnightblue        purple 
          224            99           301 

Cutoff used: 0.8 
Number of modules identified: 13 

Calculating module-module similarity based
  on module-eigengene-expression...Done.
Tidying module names...Done.
Plotting adjacency matrix for module-module similarity...

---------------------------------------------------
6. Tidy modules (clusters) 
---------------------------------------------------

# mods <- timecourseRnaseq::make_modules(
#   data = tidydata.db[[which.sample]],
#   log2 = TRUE,
#   id_column = "gene_name",
#   min_expression = NULL,  # automatically estimated
#   min_timepoints = NULL,  # automatically estimated
#   method = "wgcna",
#   qc = TRUE,
#   sim_method = "kendall",
#   soft_power = 15,      # automatically estimated
#   min_module_size = 50,
#   max_modules = 16,
#   merge_cutoff_similarity = 0.9,
#   plot_network_min_edge = 0.5,
#   plot_network = TRUE,
#   tidy_modules = TRUE
#   )

Modify network plot

Internal function; use ::: to call

trash <- purrr::map(
  sample.names,
  ~ timecourseRnaseq:::plot_adj_as_network(
      matrix = mods[[.x]][["adj_matrix_ME"]][["ME"]],
      # layout = igraph::layout.sugiyama,
      layout = igraph::layout_in_circle, # changed 
      min_edge = 0.6,
      node_label_size = 1.2,
      node_size = 20,
      edge_size = 3,
      node_frame_col = "grey20",
      node_fill_col = "grey80",
      vertex.frame.width = 3
    )
)

Visualizing a simplified representation of the circadian GCN

Visualizing a simplified representation of the circadian GCN

Visualizing a simplified representation of the circadian GCN

Annotate the network

Obtain a list of genes in each module

module_genes <- purrr::map(
  sample.names,
  ~ mods[[.x]][["module_genes"]]
) |> 
  setNames(sample.names)

Identify rhythmic modules

db_rhy <- purrr::map(
  sample.names,
  ~ load_rhy_genes(
      sample = .x
    )
) |> 
  setNames(sample.names)

###-###-###-###-###-###-###-###-
# Set your p-value of choice
# col_pval = "BH.Q"
col_pval = "default.pvalue"
# col_pval = "raw.pvalue"
###-###-###-###-###-###-###-###-
l_module_genes <- purrr::map(
  sample.names,
  ~ module_genes[[.x]] |> 
    arrange(module_identity) |> 
    group_split(module_identity) |> 
    purrr::map(
      ~ .x |> pull(gene_name)
    ) |> 
    setNames(unique(module_genes[[.x]]$module_identity))
) |> 
  setNames(sample.names)

l_rhy_genes <- purrr::map(
  sample.names,
  ~ db_rhy[[.x]] |> 
    purrr::map(
      ~ .x |> 
        filter(
          if_all(
            all_of(col_pval),
            ~ .x < 0.05
          )
        ) |> 
        filter(
          ID %in% unlist(l_module_genes)
        ) |> 
        pull(1) |> 
        unique()
    ) |> 
    purrr::compact()
) |> 
  setNames(sample.names)

Modules vs. Rhythmic genes

trash <- purrr::map(
  sample.names,
  function(x) {
    writeLines("
    #####################################################
    How many genes are in each of my geneset of interest?
    #####################################################")
    
    ## MAKE YOUR LIST OF GENES OF INTEREST ##
    
    # LIST ONE - WGCNA modules
    list1 <- l_module_genes[[x]]
    sapply(list1, length) |> print()
    
    ## LIST TWO - rhythmic genes
    list2 <- l_rhy_genes[[x]]
    sapply(list2, length) |> print()
    
    ## CHECK FOR OVERLAP
    # define size of genome
    size = length(unique(c(unlist(list1), unlist(list2))))
    # make a GOM object
    gom.1v2 <- GeneOverlap::newGOM(
      list2, 
      list1,
      genome.size = size
    )
    GeneOverlap::drawHeatmap(
      gom.1v2,
      adj.p = TRUE,
      cutoff=0.05,
      what="odds.ratio",
      # what="Jaccard",
      log.scale = T,
      note.col = "black",
      grid.col = "Oranges"
    )
  }
)

    #####################################################
    How many genes are in each of my geneset of interest?
    #####################################################
 C1  C2  C3  C4  C5  C6  C7  C8  C9 C10 C11 C12 C13 
843 257 163  78 251 347 342  82 234 173  66 174  90 
      ARS    empJTK GeneCycle       JTK    meta2d      RAIN 
      141       405       247       110       156       363 

    #####################################################
    How many genes are in each of my geneset of interest?
    #####################################################
 C1  C2  C3  C4  C5  C6  C7  C8  C9 C10 C11 C12 C13 C14 C15 
634 240  65 118 857 276 180 157 105 163  71  60 403 148 104 
      ARS    empJTK GeneCycle       JTK    meta2d      RAIN 
      468       594       293       229       484       683 

    #####################################################
    How many genes are in each of my geneset of interest?
    #####################################################
 C1  C2  C3  C4  C5  C6  C7  C8  C9 C10 C11 C12 C13 
149 224 164 373 301 159  53 424 542 413 126  99 434 
      ARS GeneCycle    meta2d 
      252       241        75